#include "stdafx.h"

namespace cpputils
{

void randomizer::seed(int Val)
{
}

int c_randomizer::max_rand() const
{
	return RAND_MAX;
}

void c_randomizer::seed(int Val)
{
	srand(Val);
}

int c_randomizer::operator()()
{
	return rand();
}

const int knuth_randomizer::MSEED = 161803398;
const int knuth_randomizer::MBIG = numeric_limits<int>::max() + 1; // overflow is not a problem

int knuth_randomizer::max_rand() const
{
	return numeric_limits<int>::max();
}

knuth_randomizer::knuth_randomizer()
{
	seed(0);
}

void knuth_randomizer::seed(int Val)
{
	int mj = MSEED - Val;
	if (mj < 0) mj += MBIG;
	ma[55] = mj;
	int mk = 1;
	for (int i = 1; i <= 54; i++)
	{
		int ii = (21 * i) % 55;
		ma[ii] = mk;
		mk = mj - mk;
		if (mk < 0) mk += MBIG;
		mj = ma[ii];
	}
	for (int k = 1; k <= 4; k++)
		for (int i = 1; i <= 55; i++)
		{
			ma[i] -= ma[1 + (i + 30) % 55];
			if (ma[i] < 0) ma[i] += MBIG;
		}
	inext = 0;
	inextp = 31;
}

int knuth_randomizer::operator()()
{
	if (++inext == 56) inext = 1;
	if (++inextp == 56) inextp = 1;
	int mj = ma[inext] - ma[inextp];
	if (mj < 0) mj += MBIG;
	ma[inext] = mj;
	return mj;
}

c_randomizer& crand()
{
	static c_randomizer Rand;
	return Rand;
}
std_randomizer& stdrand()
{
	struct init_randomizer
	{
		operator std_randomizer() const
		{
			std_randomizer Rand;
			Rand.seed(rand());
			return Rand;
		}
	};
	static std_randomizer Rand = init_randomizer();
	return Rand;
}

int uniform(int N, randomizer& Rand)
{
	long long RandBound = (long long)Rand.max_rand() + 1;
	assert(N > 0);
	long long R = Rand();
	long long Bound = RandBound;
	for (;;)
	{
		if (Bound >= N)
		{
			long long Temp = R - R % N;
			R -= Temp;
			Bound -= Temp;
			if (Bound >= N) return (int)R;
		}
		R = R * RandBound + Rand();
		Bound *= RandBound;
	}
}

double uniform(randomizer& Rand)
{
	return (Rand() + 1.0) / (Rand.max_rand() + 2.0);
}

double gaussian(randomizer& Rand)
{
	double X, Y, W;
	for (;;)
	{
		X = uniform(Rand) * 2 - 1;
		Y = uniform(Rand) * 2 - 1;
		W = X * X + Y * Y;
		if (W > 0 && W <= 1)
			return X * sqrt(-log(W) * 2 / W);
	}
}

void histogram_randomizer::set_histogram(subarray<const double> Probs)
{
	PSums.resize(Probs.size() + 1);
	PSums[0] = 0;
	for (int I = 1; I <= Probs.size(); I++) PSums[I] = PSums[I - 1] + Probs[I - 1];
}

subarray<const double> histogram_randomizer::get_partial_sums() const
{
	return PSums;
}

double histogram_randomizer::get_partial_sum(int Min, int Max) const
{
	return PSums[Max + 1] - PSums[Min];
}

int histogram_randomizer::max_rand() const
{
	return PSums.size() - 2;
}

int histogram_randomizer::operator()()
{
	return operator()(0, PSums.size() - 2);
}

int histogram_randomizer::operator()(int Min, int Max)
{
	return upper_bound(PSums.begin() + 1, PSums.end() - 1, PSums[Min] + uniform() * get_partial_sum(Min, Max)) - PSums.begin() - 1;
}

}
